function [diff, Lnew] = iterate_dist(M, L, G, Q, p_exit, yf, param,glob,options)
           
    Lnew    = Q'*((1-param.gamma).*(1-p_exit).*L) + M*G;  

    q = yf/param.C;
    diff = (sum(Lnew.*upsilonp(q,param).*q))^(-1) - param.D; 
    
end